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Abstract 

The general solution of the one-dimensional stationary Schrodinger 
equation in the form of a formal power series is considered. Its effi- 
ciency for numerical analysis of initial value and boundary value prob- 
lems is discussed. 

Consider the equation 

(pu'Y + qu = uj 2 u (1) 

where we suppose that p, q and u are complex-valued functions of an inde- 
pendent real variable x G [0,a], uo is an arbitrary complex number, p and 
q are supposed to be such that there exists a solution g of the equation 
(P9q)' + Q9o = on (0, a) such that g 6 C 2 (0,a) together with l/g are 
bounded on [0,a] and p e C 1 (0, a) is a bounded nonvanishing function on 
[0, a]. Denote g = y/pgo- In [3 J with the aid of pseudoanalytic function theory 
[2] the following result was obtained. 

Theorem 1 The general solution of (Qp has the form 

u = c\Ui + c 2 u 2 (2) 
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where c\ and c 2 are arbitrary complex constants, U\ and u 2 are defined as 
follows 



00 n 00 r 

UJ , sr^ OJ 



Ui = go Y\ -r^W and u 2 = g V —X™ (3) 



n\ * — ' n\ 

even n=0 odd n=l 



where X^ and X^ are introduced by the following recurrent equalities 

XC°) = 1, X® = 1, (4) 

and for nGN, 



x = < 



n 



n 



* (n-1) (f)0o(£K for an odd rz 



-^ <n ^(0# 2 (0^ for an even n 



(5) 



n 



?? 



^ (n_1) (0^" 2 (0^ for an odd n 



l \Cjgl(C)di for an even n 



(6) 



Another representation of the general solution of ([I]) as a formal power se- 
ries has been known since quite long ago (see [H Theorem 1]) and used for 
studying qualitative properties of solutions. The parameter uj participated in 
that representation in a very complicated manner which made that form of a 
general solution too difficult to be used for quantitative analysis of spectral 
and boundary value problems. To the contrast, the solution (J2J), Q is a 
power series with respect to uj which makes it really attractive for numerical 
solution of spectral, initial value and boundary value problems. 

The required for ([3]) particular solution g can be constructed in a similar 
way. By analogy with theorem [T] the following result can be obtained. 

Theorem 2 Suppose that q G C[0, a]. The general solution of the equation 

d 2 u(x) 



dx 2 



+ q{xju(x) 







(7) 
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on (0, a) has the form 

u = C]Ux + c 2 u 2 (8) 
where c\ and c 2 are arbitrary constants; u\, u 2 are defined as follows 



= E 



and u 2 = 



even n=0 odd n=l 

and X( n \ are introduced by the following recurrent equalities 
and /or fiGN, 



(9) 



(10) 



n 



X( n for an even 



o 



I (n) (x) = < 



n y X (n_1) (0?(0^ for an odd n 
o 

/*<-'>(fl 9 «K for an even „ 



J X^-^^dt, for an odd n 



(11) 
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Proof. First of all let us verify that both series are uniformly convergent on 
the considered interval. For this purpose we notice that for an even n, 



x: 



\ n/2 / \ n/2 

< ( max \q(x)\ ) x n < I max \q(x)\ I a n . 

- \xe[0,x] J ~ \xe[0,a] J 



Thus, the members of the series in ui can be estimated by constants: 



\x( n ){x)\ 



< 



(max^gjo ]|q(x)|) n a n { . n 

^ 1 — '- tor any x G [0, aj and the series ^ converges where c = 

1 /9 

(maXj.g^aj |o(a;)|) a. Then by the Weierstrass theorem the series in u± is 
uniformly convergent. The uniform convergence of the series in u 2 (as well 
as of the series of derivatives) can be shown by analogy. 
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Now we can apply the operator to each of the series. Note that for 
an even n > 0, 

f-X^ = n^X^ = (n - l)nqX (n ' 2 \ 
dx 2 dx 



Thus, 



d 2 ^ x("- 2 ) ^ IW 

even n=2 v 7 even n=0 



For an odd n > 1 (obviously, ^-j^r- = 0) we have 

-^iW = n ± X ^ = (n- l)nqX^ n ' 2 \ 
And consequently 

d^ U2 = q 2. = 9 2. — = ^ 2 - 

odd n=3 odd n=l 

Thus, u\ and «2 are solutions of (j7]). Last step is to verify that their 
Wronskian is different from zero at least at one point. It is easy to see that 
the Wronskian has the form 

~ X(n) \ / - X (m) \ _ / ~ X(n) \ / - ^(m) ' 

2^ n! / I 2^ m! / ^ n! ^ m! 

veven n=0 / \even m=0 / \odd n=l / \odd m=l 

At the point zero all and vanish except for X<® and Thus 
the Wronskian is equal to 1 at x — and hence the functions Ui, u<i are 
linearly independent that finishes the proof. ■ 

This theorem in different, quite more difficult notations was known al- 
ready at the beginning of the last century (see [8]) and was also used for 
qualitative analysis of solutions of 0. Nevertheless here we want to empha- 
size the extreme usefulness of theorem [2] as well as of theorem [1] for numerical 
solution of initial value and boundary value problems for second-order lin- 
ear differential equations, which escaped attention of researchers working in 
numerical solution of ordinary differential equations. First of all, let us no- 
tice that in (jHJ) ci = u(0) and C2 = u'(0). That is the representation of a 
general solution (jSJ), (Q is very convenient for solving initial value problems. 
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We remind that solution of boundary value problems for (|7j) reduces to solv- 
ing a couple of initial value problems (see, e.g., [5J), so this property of the 
representation (jSj), (Q is well suited for solving boundary value problems as 
well. 

Moreover, note that very often, e.g., in electromagnetic theory (see, [7J) 
it is necessary to solve the equation 

- — 7V + w 2 q(x)u{x) = (13) 
ax z 

for different values of the complex constant u 2 . According to theorem [2] its 
general solution can be represented as follows 



u = c\ — —. — h c 2 — r 

even n=0 odd n=l 



with and defined by (flQl)-(lT2l. Thus, once and X^ up to 
a certain order N are calculated, an approximate solution of ( fT3l) is just a 
polynomial in to with calculated coefficients X^ n ^ and X^ n \ This observation 
is also valid in the case of the solution ([2]), ([3]) of equation ([1]). This property 
is very useful for numerical solution of corresponding spectral problems which 
then reduces to finding zeros of polynomials with respect to u. In the present 
work we are more interested in studying the convergence and accuracy of the 
numerical method based on representations of the form (jSJ), (JHD in comparison 
with known standard algorithms. 

An important feature of the representation (jSJ), (Q is that it is well suited 
for symbolic calculations in principle in a general case. The coefficient q can 
be interpolated arbitrarily accurately by means of a polynomial or splines and 
then all integrations in ffTTl) and (fT2"|) can be done symbolically in a package 
like Mathematica or Maple. In the present work we made use of Matlab 7 and 
compared our results with standard Matlab ODE solvers [I], [5], especially 
with ode45 which in the considered examples gave always better results than 
other similar programs. 

Consider the following initial value problem for ([7]): q = — c 2 , w(0) = 1, 
u'(0) = — 1 on the interval (0, 1). For c = 1 the absolute error of the result 
calculated by ode45 (with an optimal tolerance chosen) was of order 10 -9 
and the relative error was of order 10~ 6 meanwhile the absolute error of the 
result calculated with the aid of theorem [2] with N (the number of formal 



5 



powers in the truncated series ([!])) from 55 to 58 was of order 10~ 16 and the 
relative error was of order 10~ 14 . Taking c = 10 under the same conditions 
the absolute and the relative errors of ode45 were of order 10~ 6 and 10 -5 
respectively meanwhile our algorithm based on theorem [2] gave values of 
order 10~ 12 in both cases. 

For the initial value problem: q = c 2 , u(0) = 1, u'(0) = —1 on the interval 
(0, 1) in the case c = 1 the absolute and the relative errors of ode45 were of 
order 10~ 8 meanwhile in our method this value was of order 10 -15 already 
for N = 50. For c = 10 the absolute and the relative errors of ode45 were of 
order 10~ 3 and 10~ 7 respectively and in the case of our method these values 
were of order 10 -11 and 10 -14 for iV = 50. 

Consider another example. Let q(x) = c 2 x 2 + c. In this case the general 
solution of ([7]) has the form 



Take the same initial conditions as before, u(0) = 1, u'(0) = — 1. Then 
meanwhile for c = 1 the absolute and the relative error of ode45 was of order 
10~ 8 and for c = 30 the absolute error was 0.28 and the relative error was of 
order 10 -6 , our algorithm (N = 58) gave the absolute and relative errors of 
order 10~ 15 for c = 1 and the absolute and relative errors of order 10 -9 and 
1Q-15 respectively for c = 30. All calculations were performed on a usual PC 
with the aid of Matlab 7. 

The results of our numerical experiments show that in fact theorem [1] 
and theorem [2] offer a new powerful method for numerical solution of initial 
value and boundary value problems for linear ordinary differential second- 
order equations. Numerical calculation of integrals involved in (jSJ), and 
in ( TIT]) . (Tl2|) does not represent any considerable difficulty and can be done 
with a remarkable accuracy. 
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